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ABSTRACT 

Many  applications  of  signal  processing  require  that  the  location  of  an  electromagnetic  emitter 
be  determined.  This  thesis  applies  an  algorithm  developed  by  Dr.  William  A.  Gardner  in  1987, 
that  geolocates  such  an  emitter.  Special  features  of  the  algorithm  not  only  allow  the  detection 
of  cyclostationary  signals  in  strong  negative  signal  to  noise  ratio  environments,  but  also 
discriminate  against  unwanted  signals  not  of  interest  that  may  either  occupy  the  same  frequency 
band  or  have  similar  signal  characteristics  as  the  signal  of  interest.  This  is  accomplished  by 
determming  spectral  correlation  densities  for  cycle  frequencies  equal  to  various  exploitable 
features  of  the  bpsk  direct  sequence  spread  spectrum  signal. 
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I.  INTRODUCTION 


A.  BACKGROUND 

Many  applications  of  signal  processing  require  that  the  direction  or  the  location  of 
the  source  of  a  signal  of  interest  be  determined.  This  requirement  could  be  for  a  variety 
of  reasons.  Active  systems  such  as  radar,  transmit  a  signal  and  evaluate  the  reflected 
energy  to  determine  first  if  a  contact  exists,  and  then  the  contact’s  bearing  and  range. 
From  this  information  contact  heading  and  speed  can  be  elicited,  and  the  contact  tracked. 
Passive  systems,  on  the  other  hand,  do  not  transmit,  but  rather  receive  signals  generated 
entirely  from  other  sources.  Submarines,  for  example,  use  passive  sonar  techniques  to 
process  noise  levels  from  all  directions  to  determine  their  exact  bearing,  and  then, 
through  various  maneuvers  of  the  ship,  an  exact  position  of  a  particular  signal  of  interest. 
Ascertaining  the  location  of  various  signals  of  interest  is  necessary  for  collision 
avoidance,  primarily,  and  tracking  and  targeting,  secondly.  This  type  of  system  uses  an 
array  to  determine  the  signal’s  direction,  or  bearing  relative  to  the  ship.  An  output  is 
produced  from  each  element  of  the  array  in  sequence  as  the  signal’s  wavefront  passes 
over  it.  The  time  delay  from  the  first  element’s  output  to  the  last  element’s  output  is 
electronically  calculated  and  a  beam  produced  in  the  direction  of  the  signal. 

The  system  described  in  this  paper  is  also  a  passive  system  but  receives 
electromagnetic  energy  rather  than  sound  pressure  energy  as  in  sonar.  And  like  the  sonar 
system  characterized,  a  time  delay  is  determined  and  converted  to  a  Line  of  Position. 
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Obtaining  two  or  more  Lines  of  Position  pinpoints  the  origin  of  the  signal  of  interest. 
This  is  called  Geolocation  when  the  source  of  the  signal  of  interest  is  on  the  Earth. 

Applications  for  this  type  of  system  are  many.  Locating  a  distress  beacon  on  the 
ground  from  a  downed  commercial  airplane  in  rugged  territory,  rescuing  a  pilot  shot 
down  over  enemy  territory,  finding  violators  of  communication  regulations, 
communication  surveillance  of  illegal  activities  such  as  drug  smuggling  operations  or  of 
hostile  forces  during  armed  conflict  or  even  during  peacetime,  and  of  course,  military 
signal  intelligence  gathering  [Ref.  3]. 

B.  OVERVIEW 

Figure  1-1  gives  a  general  idea  of  the  problem  addressed  in  this  paper.  An  emitter 
transmits  an  electromagnetic  signal  that  is  received  by  a  pair  of  receivers,  the  platform 
receiver  pair.  The  emitter  could  be  either 


a  previously  unknown  signal  source,  or  a 

known  source  that  has  changed  position. 

— 

In  either  case,  the  problem  is  first  to 

PIATPOBM  "W 
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1 
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detect  the  signal  and  then  geolocate  the 

■\  \ 

A 

source  of  its  emission. 

EMTITER 

Detecting  the  signal  is  not 

Figure  1-1  General  scenario  of  geolocation 
necessarily  a  simple  task.  Other  signals  pj.o5jen^ 

not  of  interest  and  noise  interfere  with 

and  mask  the  signal  of  interest,  as  shown  in  Figure  1-2.  The  signals  not  of  interest 
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Figure  1-2  Real-world  scenario  showing  interfering  noise  and  signals  not  of 
interest. 


might  occupy  the  same  frequency  band  as  the  signal  of  interest  and  might  also  be  at  a 
higher  power  level.  In  addition,  noise  further  hides  the  signal  of  interest  by  "cluttering" 
the  electromagnetic  spectrum,  making  it  more  difficult  to  detect  and  extract  the  desired 
signal.  Since  the  type  of  signals  that  are  of  interest  in  this  paper  are  Direct  Sequence 
Spread  Spectrum  (DSSS)  signals,  and  are  usually  imbedded  in  noise,  special  techniques 
must  be  employed  to  recover  these  signals  from  the  noise. 

In  the  following  chapters  specific  problems  of  detecting  DSSS  signals  are  discussed 
and  special  methods  that  are  used  to  overcome  them  introduced.  Then,  having 
established  the  groundwork  for  detecting  DSSS  signals,  a  method  to  geolocate  the  signal 
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of  interest,  the  Spectral  Coherence  Alignment  (SPECCOA)  algorithm,  is  described  that 
is  especially  tolerant  to  noise  and  interfering  signals  not  of  interest.  Results  of  testing 
the  SPECCOA  algorithm  using  computer  generated  test  signals  are  presented.  Varying 
degrees  of  noise  are  used  ranging  from  positive  signal  to  noise  ratios  (SNR)  to  establish 
a  baseline  reference,  to  negative  signal  to  noise  ratios,  which  means  that  the  noise  level 
is  greater  than  the  signal.  In  addition,  interfering  signals  are  inserted  to  test  the 
algorithm’s  tolerance  to  signals  not  of  interest. 
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n.  SIGNAL  CHARACTERISTICS 


This  chapter  gives  a  brief  description  of  the  characteristics  of  the  type  of  signal 
studied  in  this  paper,  a  binary -phase  shift  key  (BPSK),  direct  sequence  spread  spectrum 
(DSSS)  signal.  Some  knowledge  of  the  type  of  signal  being  studied  and  its 
characteristics  is  necessary  to  facilitate  a  better  understanding  of  the  information 
presented  in  later  chapters. 

A.  SPREAD  SPECTRUM 

Spread  Spectrum  (SS)  is  defined  as  a  signal  whose  transmitted  bandwidth  is 
substantially  wider  than  the  bandwidth  of  the  information  itself.  In  addition,  the 
increased  bandwidth  must  be  caused  by  a  spreading  signal  or  code  called  the  modulation 
which  must  also  be  known  by  the  recipient  in  order  to  despread  the  signal.  SS 
technology  has  been  in  existence  for  over  thirty  years.  A  SS  signal  has  many  features 
that  make  its  use  desirable  for  communications  in  general  and  military  applications  in 
particular.  It  provides  jamming  resistance,  interference  rejection,  multiple-access 
capability,  and  low  probability  of  intercept  (LPI)  capability.  [Ref  6;  p.  404][Ref  8:  p. 
2-1] 

Using  SS  techniques,  the  energy  from  a  signal  is  spread  over  a  wide  band  of 
frequencies  by  means  of  a  spreading  code.  The  spreading  before  transmission  and 
subsequent  despreading  at  the  receiver  provides  a  spread  spectrum  signal  with  its 
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desirable  characteristics.  Orthogonal  spreading  codes  allow  numerous  signals  to  be 
transmitted  simultaneously  on  the  same  frequency  band.  A  potential  receiver  selects  the 
desired  signal  by  despreading  with  the  correct  code  while  other  signals,  because  of  their 
unique  codes,  are  left  undisturbed.  This  is  code  division  multiple  access  (CDMA). 
Interference  with  narrow-band  signals  is  minimized  because  the  SS  signal’s  power  is 
spread  over  such  a  large  bandwidth  that  the  in-band  narrow-band  signal  to  spread 
spectrum  signal  ratio  is  very  large.  This  explains  the  other  desirable  features  of  SS 
signals.  Since  the  power  level  is  very  low  at  each  frequency  over  which  the  signal  is 
spread,  the  signal  has  a  low  probability  of  detection  and/or  intercept  (LPD/LPI),  and  an 
interfering  or  potential  jamming  signal  is  spread  during  the  despreading  process  of 
recovering  the  desired  signal.  This  gives  the  signal  of  interest  a  high  processing  gain 
which  necessitates  a  much  stronger  interfering  or  potential  jamming  signal  in  order  to 
prevent  reception  of  the  signal  of  interest. 

Some  categories  of  Spread  Spectrum  are: 

•  Frequency  Hop  (FH) 

•  Stacked  Carrier 

•  Direct  Sequence  (DS) 

•  Hybrid 

Frequency  Hop  Spread  Spectrum  transmits  communication  signals  by  rapidly 
hopping  between  several  carrier  frequencies,  numbering  from  just  a  few  to  many 
thousands.  These  systems  can  be  either  a  fast  or  a  slow  hopper,  depending  on  whether 
several  data  bits  are  transmitted  per  hop  frequency  (slow  hopper)  or  several  hop 
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frequencies  are  used  per  data  bit  (fast  hopper).  This  method  makes  the  signal  more  jam 
resistant.  In  addition,  it  is  difficult  to  collect  the  entire  signal,  giving  it  a  low  probability 
of  intercept.  [Ref.  8] 

Stacked  carrier  signals  provide  extremely  reliable  communications  by  transmitting 
simultaneously  on  two  or  more  frequencies.  This  increases  anti-jamming  resistance  and 
lowers  the  error  rate  of  transmission.  Two  drawbacks  of  using  stacked  carrier  are  its 
high  signal  to  noise  ratio  making  it  susceptible  to  interception  and  its  inefficient  use  of 
the  frequency  spectrum. 

Direct  Sequence  Spread  Spectrum  uses  a  pseudo-random  spreading  sequence  at  a 
chip  frequency,  f^,  to  modulate  the  carrier 
(frequency  fj,  and  spread  a  signal  over  a 
wide  frequency  band.  See  Figure  2-1. 

The  degree  of  spreading  is  a  function  of 
the  chip  frequency  and  can  be  determined 
by  dividing  the  chip  frequency  by  one- 
half  of  the  original  bandwidth,  or 
spreading  factor  =  2f(./BWo,  where  BW<, 
is  the  bandwidth  of  the  original  signal. 

Note  also  that  the  power  is  spread  over  a  wide  frequency  band  and  the  peak  power  is 
lowered,  potentially  to  below  that  of  the  noise  threshold.  Thus,  the  spreading  can  be  to 
such  a  degree  that  any  background  noise  masks  the  final  transmitted  signal.  Using  this 
method  makes  detection  of  the  signal  more  difficult  because  an  ordinary  receiver 


Figure  2-1  Illustration  of  chip  features 


spreading  original  signal  and  lowering  its 
peak  power. 
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perceives  only  noise.  In  addition,  the  signal  is  harder  to  jam  since  a  Jammer  would  have 
to  Jam  numerous  frequencies  and  the  process  of  despreading  the  signal  tends  to  spread 
the  Jamming  signal  as  described  previously. 

The  Hybrid  spread  spectrum  signal  uses  a  combination  of  other  spread  sequence 
techniques,  the  most  common  of  which  is  the  frequency  hop/direct  sequence  hybrid. 
Frequency  hopping  is  used  to  achie\'e  the  same  large  bandwidth  of  a  direct  sequence 
signal  but  with  a  lower  chip  rate.  This  is  desirable  since  it  is  more  difficult  to  achieve 
the  large  bandwidth  using  direct  sequence  methods  alone.  [Ref.  8] 

B.  BINARY  PHASE  SHIFT  KEYING  (BPSK) 

A  BPSK  signal  is  produced  by  shifting  the  phase  of  a  sinusoidal  carrier  180°  with 
a  binary  data  sequence.  See  Figure  2-2.  Shown  is  the  sinusoidal  carrier  signal,  cos(a>,t), 
and  the  unipolar  data  sequence,  d(t),  which  is  shifted  to  a  bipolar  signal,  m(t).  to 
modulate  the  carrier,  resulting  in  the  signal  to  be  transmitted,  m(t)cos(ajet).  The  effect 
of  the  modulating  signal  is  to  shift  the  phase  of  the  carrier  by  1 80^  each  time  the  data  bit 
changes  from  a  binary  1  to  a  binary  0,  or  vice  versa.  [Ref.  6;  pp.  332-336]  The 
number  of  data  bits  per  second  is  the  data  rate  or  bit  frequency,  f^,  of  the  signal.  Also 
of  note  is  the  carrier  frequency  of  the  signal,  f^  i=(j}J2ir).  The  bandwidth  of  the 
transmitted  signal  is  then  [(fo+fb)  -  (fo-fb)]  or  2fb.  In  the  next  section  it  is  shown  how 
BPSK  and  DSSS  are  combined  to  produce  a  signal  with  a  substantially  increased 
bandwidth. 
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carrier  signal,  cos(cOet) 


modulated  carrier  signal,  m(t)cos(cOct) 
Figure  2-2  BPSK  signal  production. 


C.  BPSK-DSSS  SIGNAL 


Important  features  of  a  BPSK-DSSS  signal  are  the  carrier  frequency,  f,^,  the  chip 
rate  feature  or  spreading  code  rate,  f ,  and  the  bandwidth  of  the  transmitted  signal.  A 
spreading  code  similar  in  appearance  to  that  of  the  binary  data  sequence  shown  in  Figure 
2-2,  is  modulo-2  added  to  the  data  sequence  and  shifted  to  obtain  a  bipolar  binary 
modulating  signal  for  the  carrier  frequency.'  Differences  between  the  data  sequence 
and  the  spreading  code  are  the  frequency  (f^  is  substantially  greater  than  f^)  and  the 
randomness  of  the  binary  sequence  -  the  best  spreading  codes  are  those  that  are  as  near 
random  as  possible,  and  of  course  the  data  sequence  is  not  random  since  it  contains  the 
information  to  be  transmitted. 

The  result  is  a  transmitted  signal  with  the  same  carrier  frequency  as  before,  f,,,  but 
with  a  significantly  increased  bandwidth  equal  to  [(fo+fc)  -  (fo-fc)]  or  2f,.  The  ratio  of 
the  new  bandwidth  to  the  original  bandwidth  is  f^/f,,  which  is  called  the  processing  gain. 
There  is  also  a  reduction  of  the  peak  power  of  the  BPSK-DSSS  signal  as  compared  to  the 
original  signal  (see  Figure  2-1,  above)  by  the  same  factor.  It  is  called  the  processing 
gain  since  it  is  the  increase  in  power  that  the  signal  gains  during  the  recovery  process. 

D.  SUMMARY 

This  chapter  describes  the  various  features  associated  with  the  signal  of  interest 
addressed  in  this  paper.  Special  intercept  techniques  are  needed  to  detect  signals  such 


'Modulo-2  is  XOR  addition;  0©0  =  0,  001  =  I.  100  =  1,  101  =0. 
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as  these  that  are  below  the  noise  threshold,  W.  A.  Gardner  in  a  good  survey  article 
[Ref.  3]  discusses  the  exploitation  of  spectral  redundancy  in  cyclostationary  signals  that 
allow  processing  signals  that  are  deeply  imbedded  in  noise.  This  property,  which  is 
inherent  to  manmade  communication  signals,  enables  a  sifting  through  of  extraneous 
noise  and  interfering  signals  to  piece  together  and  recover  the  desired  signal.  [Ref.  4] 
introduces  several  methods  that  use  cyclostationarity  to  geolocate  a  signal  of  interest. 
These  items  are  discussed  in  more  detail  in  Chapter  3. 
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ffl.  SPECTRAL  CORRELATION 


Cyclostationarity  is  a  spectral  redundancy  feature  that  allows  the  detection  of 
signals  deeply  imbedded  in  noise,  especially  those  signals  that  are  Direct  Sequence 
Spread  Spectrum  signals.  This  chapter  covers  the  theory  associated  with  this  property 
and  how  it  is  used  to  uncover  hidden  signals. 

A.  DETECTING  HIDDEN  PERIODICITY 

This  theoretical  development  is  taken  from  W.  A.  Gardner’s  Statistical  Spectral 
Analysis,  A  Nonprobabilistic  Theory,  [Ref.  7,  pp.  359-362],  which  can  be  studied  for  a 
more  detailed  treatment  of  the  subject.  Spectral  lines  can  be  generated  by  putting  a 
signal  through  a  quadratic  time-invariant  (QTI)  transformation,  e.g.,  H[s(t)s(t  -  r)]. 

Starting  with  a  signal  x(t)  which  contains  a  finite-strength  additive  sinewave 
component  (an  ac  component)  with  frequency  a 

a  cos(2Ttaf  +  0),  0, 


if  the  Fourier  coefficient 


M 


lim 

T  -  c 


1  ^ 

-  f  dt 


(3.2) 


-m 


exists  and  is  not  zero,  then 
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(3.3) 


m;  =  -ae‘^. 

2 

In  this  instance,  the  power  spectral  density  of  x(t)  includes  a  spectral  line  at  frequency 
f  =  a  and  its  image  f  =  -a.  In  other  words  the  power  spectral  density  contains  the 
additive  term 

|M;!^  m-  a)  -  6(f  ^  a)]  (3.4) 


where  6(*)  is  the  impulse  function.  See  Figure  3-1.  The  spectral  lines  of  x(t)  at  f  -  ±a 
are  due  to  the  finite  additive  sine  wave  component.  [Ref.  8][Ref.  3,  pp.  16.17] 

Rewriting  x(t)  in  terms  of  (3. 1)  and  lumping  all  else  that  x(t)  may  contain  into  n(t) 

gives 

x(t)  =  a  cos(27caf  ♦  0)  +  n(t)  (3*5) 


where  n(t)  is  random.  The  signal  x(t)  can 
be  said  to  have  hidden  periodicity  if  n(t) 
is  much  greater  than  the  sine-wave 
component  making  it  not  readily 
discernable  with  a  perfunctory 
observation.  However,  due  to  the 
spectral  lines  at  f  =  ±a  the  hidden 

Figure  3-1  Graph  of  the  power  spectral 
periodicity  can  be  detected  by  examining  density  with  spectral  lines  at  ±  a. 

the  signal  more  carefully.  [Ref.  3,  pp.  16.  17] 
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This  signal  has  first-order  periodicity  with  frequency  a  and  its  spectral  components 
can  be  detected  without  any  special  transformation  of  the  original  signal.  Many  man¬ 
made  signals,  however,  contain  hidden,  higher-order  periodicities  that  do  not  generate 
spectral  lines  at  the  frequencies  of  the  hidden  periodicities.  Using  a  simple  squarer 
quadratic  transformation,  i.e,  x-(t)  =  x(t)x(t),  some  signals  can  be  changed  into  a  first- 
order  periodicity  exposing  their  hidden  spectral  components.  For  other  signals,  the 
simple  squarer  transformation  does  not  work  and  a  better  transformation  is  one  that  uses 
a  time  delay  such  as 

y(f)  =  jc(r)jc(r  -  t) 

for  various  nonzero  delays  t  [Ref.  3,  p.  18],  This  delayed  transformation  gives  rise  to 
the  cyclic  autocorrelation  fiinction.  It  can  be  shown  [Ref.  7,  chap.  10]  that  the  signal 
x(t)  contains  second-order  periodicity  at  frequency  a  if 


exists  and  is  not  zero  for  some  non-zero  a.  is  defined  as  the  cyclic  autocorrelation 

of  x(t)  for  cycle  frequency  parameter  a.  The  definition  of  this  function  comes  from  the 
Fourier  coefficient  of  the  additive  sinewave  component  with  cycle  frequency  a. 
Substituting  (3.6)  into  (3.2)  gives 

My  =  <y(t)e-^’“‘>  =  <x(t)x(t  -  t)e72<t«»>^  (3.8) 
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where  <  •  >  is  the  time  average 


lim 

r  -  oo 


(3.9) 


Adjusting  the  arguments  for  symmetry  yields 

A/;  =  <y(r)c  -^2««t>  ^  ^  (3,10) 


which  leads  to  the  definition  of  /?“(t)  shown  in  (3.7).  Taking  the  cyclic  autocorrelation 

function  one  step  further  leads  to  the  spectral  correlation  density  and  the  development  of 
the  bifrequency  plane. 

B.  BIFREQUENCY  PLANE 

From  statistical  processing  it  is  known  that  the  power  spectral  density  and  the 
autocorrelation  function  are  fourier  transform  pairs,  i.e., 

00 


and 


K,(t)  " 


(3.12) 


Similarly,  it  can  be  shown  [Ref.  7,  p.  389]  that  the  cyclic  spectral  correlation  density  and 
the  cyclic  autocorrelation  function  are  also  Fourier  transform  pairs. 


15 


5;(/)  = 


(3.13) 


and 

<a 

Rxi^)  =  { df  (3.14) 

•>00 

where  a  is  the  cycle  frequency.  Using  this  concept  a  correlation  can  be  taken  in  the 
frequency  domain  of  a  signal  x(t).  The  correlation  parameter  is  called  the  cycle 
frequency  a,  and  a  plot  of  the  frequency  vs.  the  cycle  frequency  produces  the 
bifrequency  plane.  Figures  3-2  and  3-3  show  how  this  plane  is  developed.  Starting  with 
the  frequency  domain  signal  shown  in  Figure  3-2,  a  correlation  is  taken  between 
frequencies  (f+a/2)  and  (f-a/2)  as  a,  the  cycle  frequency,  varies  from  -f,  to  -(-fj,  where 
f,  is  the  sample  frequency.  Mathematically  the  product  X(f  +  al2)X\f  -  a/2)  is  taken 
and  the  magnitude  plotted  on  the  bifrequency  plane  shown  in  Figure  3-3,  at  coordinates 
(f,  a).  This  unsmoothed  bifrequency  plane  is  called  the  cyclic  periodogram.  The 

unsmoothed  bifrequency  plane  may  be 
smoothed  either  in  time  or  frequency  to 
minimize  variance  effects  and  noise  (Ref.  8: 
Ch.  13].  Figure  3-4  is  a  plot  of  a  typical 
bifrequency  plane  of  a  BPSK  signal  with 
carrier  frequency  f^.  Important  features  on 
this  plot  are  the  power  spectrum  shown  for 


X(l| 


Figure  3-2  Frequency  domain  signal 
showing  correlation  being  taken  at 
(f  ±  a/2). 
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cycle  frequency  equal  to  zero,  the 
spreading  code  or  chip  rate  features  (the 
four  smaller  peaks),  and  the  carrier 
feature,  the  targe  peak  occurring  at  cycle 
frequency  equal  to  twice  the  carrier 
frequency,  or  a  =  2fo. 

Figure  3-3  General  shape  of  the 
bifrequency  plane.  Magnitudes  of  the 
C.  APPLICATIONS  correlation  are  plotted  at  coordinates 

(f,  a)  as  f  varies  over  the  signal’s 
There  are  numerous  applications  of  bandwidth,  and  a  varies  from  -2fs  to  +2fs. 

spectral  correlation.  Some  of  these  [Ref.  8:  pp.  30-34]  are: 

•  Detection  and  classification  of  signals. 

•  Parameter  estimation  such  as  carrier  frequency  and  keying  rate. 

•  Spatial  filtering. 

•  Direction  fmding. 

•  Frequency  shift  filtering  for  signal  extraction. 

•  Frequency-shift  prediction. 

•  Time  difference  of  arrival. 

[Ref  8]  gives  further  information  on  specific  applications.  The  last  item,  time  difference 
of  arrival,  is  the  subject  of  this  thesis.  Using  the  spectral  correlation  densities  for  the 
cycle  frequency  not  equal  to  zero,  as  shown  in  Figure  3-4,  lines  of  position  are 
developed  and  the  location  of  the  emitter  of  the  signal  determined.  Determining  a  line 


a 
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IV.  TDOA 


A.  OBTAINING  A  FIX 

In  Naval  terms,  to  obtain  a  fix  means  to  determine  a  ship  s  position.  Several 
methods  exist  that  are  used.  One  such  system  developed  during  World  War  n.  is  the 
LORAN  electronic  navigation  system  and  is  briefly  mentioned  here  because  its  method 
of  fixing  a  ship's  position  is  similar  to  the  method  used  in  this  thesis  to  geolocate  an 
emitter.  LORAN,  from  long  range  navigation,  produces  hyperbolic  lines  of  position 
based  on  the  time-delay  difference  between  signals  received  from  a  pair  of  land-based 
transmitting  stations  [Ref  1:  pp.  345-348].  Figure  4-1  depicts  this  situation.  A  pair  of 

fixed  transmitters  separated 
by  a  known  distance  along 
the  baseline  transmits 
signals  that  produce 
hyperbolic'  isochrones. 
Isochrones  are  plotted  on 
special  navigation  charts 
and  the  correct  isochrone 
is  determined  based  on  the 


A  hyperbola  is  a  mathematical  term  defined  as  the  locus  of  points  whose  difference  in  distance  from  two  fixed 
points  is  constant. 
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transmitter  pair's  unique  identification  number  and  the  time  of  arrival  difference  between 
the  received  signals,  i.e..  the  At  or  difference  between  t,  and  t,.*  A  ship  receives  these 
signals  in  an  attempt  to  fix  its  position.  Possible  ship  positions  could  be  anywhere  along 
the  indicated  hyperbola;  three  possible  locations  are  marked.  This  establishes  a  single 
line  of  position.  Repeating  the  procedure  for  a  second  pair  of  transmitters  gives  a  second 
line  of  position.  The  intersection  of  the  two  hyperbolic  isochrone  fixes  the  ship's 
position,  as  in  Figure  4-2. 

B.  GEOLOCATION 

In  the  above  scenario,  the  location  of  the  transmitters  is  already  known  since  their 
positions  are  printed  on 
navigation  charts. 

However,  for  the  problem 
addressed  in  this  thesis,  the 
situation  is  reversed,  that 
is  to  say,  the  signal  is  a 
single  point  origin,  rather 
than  a  pair  of  transmitters, 
and  the  receivers  are  a  pair 

‘The  two  transmitters  are  designated  Master  and  Slave  since  the  Slave  transmitter  does  not  transmit  until 
the  Master  station  transmission  reaches  the  Slave  station.  This  avoids  position  ambiguities  near  the  line  exactly 
between  the  two  transmitters.  If  both  stations  transmitted  simultaneously,  there  would  be  zero  time-delay  difference 
exactly  on  the  center  line,  and  if  the  ship  was  not  exactly  on  the  center  line,  ambiguity  would  exist  as  to  which  side 
of  the  line  the  ship  was  actually  located,  and  a  position  error  would  result. 


Figure  4-2  Intersection  of  two  Isochrones  to  get  a  Fix. 


20 


of  antenna  platforms  separated  by  a  known  distance.  See  Figure  4-3.  In  this  situation, 
each  platform  receives  the  same  signal  from  an  unidentified  source  but  at  a  different 
time.  At  time  t^,,  if  the  signal  at  receiver  platform  one  is  x(t),  and  the  signal  at  receiver 
platform  two  is  y(t),  then  y(t)  =  x(t  -  d)  where  d  is  the  delay  or  difference  of  times 
between  the  signal  arriving  at  platform  one,  x(t),  and  platform  two,  y(t).  The  difference 
yields  a  hyperbola  of  possible  emitter  locations,  the  hyperbolic  line  of  position  at  time 
to.  A  second  hyperbolic  line  of  position  is  obtained  a  short  time  later  at  time  t,  as  the 
platform  receiver  pair  moves  a  short  distance  from  its  original  position.  The  intersection 
of  the  two  hyperbolas  geolocates  the  emitter  source. 

Appendix  A  provides  a  mathematical  proof  to  show  that  the  above  described 
scenario  does  indeed  produce  a  hyperbola.  Next,  Chapter  V  describes  the  algorithm  that 
produces  the  time  delay  to  obtain  the  hyperbolic  lines  of  position  just  discussed. 


Figure  4-3  Moving  antenna  platforms  Geolocating  an 


emitter. 
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V.  SPECCOA  ALGORITHM 


A,  INTRODUCTION 

In  Chapter  HI,  it  is  shown  how  a  spectral  redundancy  feature  called 
cyclostationarity  is  used  to  detect  signals  that  are  deeply  imbedded  in  noise,  i.e.,  signals 
with  a  negative  signal-to-noise  ratio  (SNR).  Once  the  signal  is  detected  by  the  platform 
receiver  pair  it  can  be  further  processed  by  Time  Difference  of  Arrival  (TDOA) 
techniques  to  geolocate  the  new  signal.  Using  the  cyclostationary  properties  of  the  signal 

ft 

of  interest,  the  auto-spectral  correlation  density  for  the  signal  at  the  first  receiver, 
5“(/),  and  the  cross-spectral  correlation  density  for  the  signals  at  both  receivers, 

5y“(/),  for  a  specific  cycle  frequency,  a,,,  are  obtained.  These  three  elements  are  the 

main  ingredients  needed  for  use  in  the  SPECtral  Coherence  Alignment  (SPECCOA) 
TDOA  algorithm  which  is  the  implementation  of  the  following  equation: 

D  =  Re  f 

-772 

Here  5  “(/)*,  denotes  the  conjugate  of  the  auto-spectral  correlation  of  the  signal  x. 


(5.1) 
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Describing  the  SPECCOA  algorithm  in  words,  the  product  of  the  two  spectral 
correlation  densities  is  inverse  Fourier  transformed.  The  result  is  multiplied  by  a 

correction  factor,  the  last  term  of  equation  (5.1).  to  ensure  the  maximum  value 

of  the  algorithm  occurs  at  the  time  delay,  d.  Finally,  the  maximum  value  of  the  real  part 
is  taken.  This  value  is  the  output  of  the  algorithm  and  is  the  delay  between  when  the 
signal  of  interest  arrives  at  receiver  one,  x(t),  and  when  it  arrives  at  receiver  two.  y(t). 
i.e.,  the  delay,  d,  in  y(t)  =  x(t  -  d).  See  Figure  5-1  below.  Depicted  is  a  typical  plot 


Figure  5-1  Typical  plot  of  SPECCOA  algorithm  output, 
of  equation  (5.1).  A  negative  lag  is  shown;  a  positive  lag  indicates  that  the  signal  first 
arrived  at  the  other  receiver.  A  lag  of  51,  meaning  51  samples  since  the  signal  is 
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digitally  processed,  is  shown  for  this  example.  The  maximum  lag  value  is  converted  to 
an  actual  time  delay  by  multiplying  the  lag  value  by  the  sampling  period.  Using  the  time 
delay,  a  hyperbola  is  generated  representing  the  line  of  position  as  described  in  Chapter 
IV.  A  second  time  delay  from  the  algorithm  generated  a  short  time  later  produces  a 
second  hyperbolic  line  of  position.  The  intersection  of  these  two  lines  of  position 
geolocates  the  emitter  of  the  signal  of  interest. 

B.  THEORETICAL  DEVELOPMENT 

This  section  shows  that  the  output  of  the  SPECCOA  algorithm  of  equation  (5.1) 
is  indeed  the  desired  time  delay,  d,  between  the  signals  x(t)  and  y(t)  =  x(t  -  d).  The 
proof  begins  by  starting  with  an  equivalent  expression  for  the  cross-spectral  correlation 

density  [Ref.  8:  ch.  7],  and  making  appropriate  substitutions. 


O/)  =  /  2)  (5.2) 

where  1/Af  is  the  size  of  the  FFT  and  indicates  frequency  smoothing.  Since  y(t)  =  x(t  - 
d),  and  using  f'  to  distinguish  from  the  argument  /  in  (5.2),  then 

'5-3) 

which  is  a  property  of  the  Fourier  transform.  Substituting  (5.3)  into  (5.2)  gives 
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(5.4) 


v.(/)  =  Wo/"  “c/2)  w./-“<y2)  ^ 


or 


-flnif*  <iJ2)d 


(5.5) 


since  by  definition,  =  X  X  ’ ,  but  here  it  is  for  a  specific  alpha,  a^. 

Using  equation  (5.5),  the  inverse  Fourier  Transform  of  the  cross- spectral  and  auto- 
spectral  correlations  is  taken. 


few  Cw} 

«0 

=  ]■  ^  (*•<«> 


-/2it>a  g/2n/c^ 


-jna^d 

e 


(5.6b) 


The  bracketed  term,  an  inverse  Fourier  transform,  has  a  maximum  value  when  the 
variable  factor,  maximum.  This  occurs  when  the  lag,  r, 

is  equal  to  the  delay,  d.  Finally,  multiplying  the  last  term  of  (5.6b)  bye-'" “^produces  the 
desired  SPECCOA  equation.  This  last  product  produces  a  second  variable  term, 
gjnaix  -d)  jj  maximum  when  the  lag  is  equal  to  the  delay.  Thus,  the  output  of 
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the  SPECCOA  algorithm  is  the  delay  between  signals  x(t)  and  y(t).  [Ref  4:  pp.  1174. 
1177] 


C.  ALGORITHM  TESTING 

1.  SSPI  Software  Package 

The  SPECCOA  algorithm  was  tested  using  simulated  signals  generated  by  a 
software  package  from  Statistical  Signal  Processing,  Inc.  (SSPI).  This  is  a 
comprehensive  package  that  generates  a  variety  of  simulated  signals.  Input  values  allow 
setting  of  sample  size,  cycle  frequency,  carrier  frequency,  desired  signal  delay,  signal 
power,  noise  power,  and  several  other  items  not  listed  here.  *  The  simulated  signals  are 
processed  by  the  software  package  to  obtain  spectral  and  cross-spectral  correlation 
densities  for  the  user  specified  cycle  frequency.  The  results  are  stored  in  data  files  for 
further  processing  at  the  users  discretion  or  for  plotting  using  the  SSPI  sxaf_plot 
command.  For  the  application  of  this  paper,  the  data  for  the  spectral  and  cross-spectral 
correlation  densities  was  input  into  the  SPECCOA  algorithm. 

2.  Testing  Overview 

Testing  was  performed  on  the  exploitable  features  shown  in  Figure  5-2;  the 
spreading  frequency  or  chip  rate  feature,  and  the  carrier  frequency  feature.  These 
features  are  exploited  by  setting  the  cycle  frequency,  equal  to  the  specific  feature  to 
be  exploited.  Thus,  to  exploit  the  chip  rate,  the  four  smaller  peaks,  is  set  to  f  or  to 


'  The  complete  list  of  items  is  given  in  Appendix  E. 
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2f„±f ,  and  to  exploit  the  carrier 
frequency,  f^,  the  larger  peak, 
is  set  to  2f„.  Since  f,  +  f  and  f^,  - 
f  are  essentially  the  same,  a  total 
of  three  sets  of  tests  were 
perfonned;  =  f,,  =  2f,„  and 

=  f,  +  f,.  Each  set  and  their 
results  is  described  next.  One 
note  that  is  pertinent  to  each  set  of 
tests  is  the  SNR  input.  Generally, 
the  SNR  referred  to  is  the  ratio,  in 
dB,  of  the  signal  input  power  to  the  total  noise  power  which  is  a  distinct  input  in 
generating  the  test  signal  as  described  in  Appendix  E.  The  final  SNR  is  actually 
somewhat  larger  since  the  broadband  noise  is  essentially  filtered  out  when  the  SPECCOA 
algorithm  calculates  a  delay,  and  only  the  inband-noise  affects  the  final  SNR.  The  final 
SNR  can  be  detennined  by  calculating  the  signal  power  to  in  band  noise  ratio:  SNR,,„,,  = 
(P,/P„)^f,  where  P,  and  P„  are  the  input  signal  and  noise  power  inserted  into  the  signal 
generation  program,  and  f  is  the  digital  spreading  code  rate  or  chip  feature.  Thus  for 
an  input  SNR  of  -15dB  and  a  f,.  of  0.0625  the  actual  in-band  SNR  is:  -0  0625 


Figure  5-2  Bifrequency  plot  of  exploitable  cycle 
frequency  features. 


or 


31.6 


(16)  =  0.506  or  -3dB  which  means  that  the  inband  SNR  is  l2dB  greater  than 


the  input  total  SNR  for  a  chip  feature  of  0.0625.  A  summary  of  the  increase  in  SNR 
over  the  input  broadband  SNR  for  all  the  chip  frequencies  used  is  shown  in  Table  5-1. 


TABLE  5-1 

SNR  INCREASE  OF  INBAND  SNR  OVER 
INPUT  TOTAL  SNR  FOR  ALL  CHIP  RATES 
USED  IN  THE  TEST  SIMULATIONS 


SNR 

INCREASE 

SNR 

INCREASE 

0.015625 

IgdB 

0.2625 

5.8dB 

0.03125 

15 

0.032 

4.9 

0.0625 

12 

0.3225 

4.91 

0.125 

9 

0.38 

4.2 

0.20 

7 

0.3825 

4.17 

0.26 

5.9 

0.4025 

4.0 

'xx/T, 


3.  Equals  Chip  Rate 

Table  5-2  shows  the  variables  input  into  the  signal  generation  program  along 
with  a  sample  of  typical  values  used.  Other  variables,  shown  in  the  complete  list  of 
Appendix  E,  were  not  changed.  Testing  consisted  of  varying  the  delay,  percent 
smoothing,  chip  feature  and  the  broadband  SNR.  The  sample  length  is  listed  in  Table 
5-2  even  though  it  was  not  a  variable  tested,  because  it  was  varied  in  response  to 
changing  the  SNR  in  order  to  obtain  an  acceptable  output.  Complete  results  of  this  test 
set  are  presented  in  Appendix  B. 
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Table  5-2  INPUT  VARIABLES  INTO  SIGNAL  GENERATION  PROGRAM 


input 

% 

no.  of 

chip  feature. 

carrier 

SNR. 

delay 

smoothing* 

freq 

samples* 

f  ^ 

‘c 

frequency, 
f  ^ 

‘0 

dB 

101 

0.01 

2048 

0.0625 

- 

0.16 

-5 

'AtAf  =  (%  smooth/ 100)(no.  of  freq  samples) 

'input  is  power  of  two;  samples  going  to  speca  program  is  (pwr  of  2  input)(  1  -  04,) 
'f /f 

I,/  1, 


One  particular  result  of  note  is  the  percent  smoothing  variation.-  In  order 
to  first  detect  a  signal  for  further  examination,  a  significant  amount  of  smoothing  is 
required,  on  the  order  of  5%.  However,  test  results  show  that  very  little  smoothing  is 
required  of  the  signal  for  insertion  of  the  spectral  correlation  density  into  the  SPECCOA 
algorithm.  Figures  5-4  through  5-7  show  the  effects  of  5 %  and  0.01  %  smoothing  (AtAf 
=  1  and  410,  respectively)  for  input  total  SNR’s  ranging  from  OdB  to  -15dB.  Clearly, 
to  detect  a  signal  and  determine  if  any  features  are  available  for  further  exploitation  a 
high  degree  of  smoothing  is  required.  Examining  the  0.01  %  smoothing  plots  would  not 
reveal  whether  anything  was  exploitable.  But.  for  input  into  the  SPECCOA  algorithm, 
the  0.01  %  smoothing  is  more  desirable.  Figure  5-8  shows  the  effect  on  the  algorithm's 
performance  for  increasing  smoothing.  Even  smoothing  as  little  as  0.39%  causes 
erroneous  results.  Inband  SNR  for  Figures  5-4  through  5-8  is  12dB  higher  than  the  input 
SNRs  shown  as  explained  earlier. 

A 

percent  smoothing  =  — 'lOO 
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General  conclusions  for  results  of  this  first  test  set  are  that  the  algorithm 
worked  well  for  varying  parameters  of  delay,  chip  feature  and  carrier  frequency  for  input 
SNRs  down  to  -15dB  corresponding  to  an  inband  SNR  of  -3dB  for  a  chip  feature  of 
0.0625/T,.  Lower  SNRs  can  be  achieved  by  increasing  the  sample  length,  see  Figure  5- 
9,  but  of  course,  the  time  to  obtain  an  output  from  the  algorithm  increases  due  to  the 
amount  of  calculations  necessary.  Thus,  the  cycle  frequency  equal  to  the  chip  rate 
feature  is  an  exploitable  feature. 

4.  Equals  Twice  the  Carrier  Frequency 

The  results  of  this  test  set  closely  resemble  those  of  set  number  one.  One 
major  difference  is  the  overall  smoother  appearance  of  the  output  graphs.  This  is  due 
to  a  greater  degree  of  zero-padding  required  to  reach  the  next  higher  power  of  two  before 
inverse  Fourier  transforming.  Complete  results  of  this  test  set  are  presented  in  Appendix 
C.  General  conclusions  for  this  test  set  are  the  same  as  for  set  number  one:  cycle 
frequency  equal  to  twice  the  carrier  frequency  is  an  exploitable  feature  for  BPSK  signals, 
but  because  of  the  extra  smoothing  due  to  zero-padding,  results  are  more  accurate. 

5.  a„  Equals  Twice  the  Carrier  Frequency  Plus  the  Chip  Rate 

The  results  of  this  test  set  are  presented  in  Appendix  D.  General  results  are 
similar  to  those  described  above,  however,  a  different  output  appearance  is  evident.  This 
is  due  to  truncation  down  to  the  next  lower  power  of  two  before  inverse  Fourier 
transforming,  and,  as  with  the  second  test  case,  the  results  are  generally  more  accurate 
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than  test  case  one.  Overall  conclusion:  cycle  frequency  equal  to  2f(,+f,  is  suitable  for 
exploitation  of  BPSK  signals. 

6.  Interference  of  Other  Signals 

The  testing  of  the  algorithm’s  ability  to  discriminate  against  signals  not  of 
interest  was  done  using  four  types  of  interference:  Multiple  signal  interference,  wide¬ 
band  signal  interference,  coband  interference,  and  narrow-band  interference.  Figure  5-3 
shows  how  the  spectral  correlation  is  able  to  distinguish  between  different  signals. 
Shown  is  the  wide-band  interferer  case.  Compare  with  Figure  5-2.  The  parameters  of 
the  interfering  signals  are  shown  in  Table  5-3,  and  Figure  5-9  shows  the  output  results. 
The  algorithm  was  able  to  sift  through  extraneous  signals  not  of  interest  and  correctly 
determine  the  delay. 
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TABLE  5-3  SIGNAL  PARAMETERS 
FOR  SIGNAL  INTERFERENCE  TEST  CASE 


.signal 

power 

carrier 

frequency' 

OdB 

0.201 

OtIB 

0.121 

OdB 

0.312 

OdB 

0.21875 

OdB 

0.16 

OdB 

0.128 

OdB 

0.16 

MULTIPLE 

INTERFERENCE 


WIDEBAND 

INTERFERENCE 


COBAND 

INTERFERENCE 


NARROW-BAND 

INTERFERENCE 


SIGNAL  OF  INTEREST 


input 

delay 


0.0196 


0.0625 


2 


Figure  5-4  Effect  of  smoothing  on  bifrequency  plane  output  and  detection  of 
exploitable  cycle  frequency  features:  a)  .01%  smoothing,  b)  5%  smoothing.  OdB 
SNR  input,  4-12dB  in-band  SNR. 


Figure  5-5  Effect  of  smoothing  on  bifrequency  plane  output  and  detection  of 
exploitable  cycle  frequency  features:  a)  .01%  smoothing,  b)  5%  smoothing.  -5dB 
SNR  input,  +7dB  in-band  SNR. 


Figure  5-6  Effect  of  smootiiing  on  bifrequency  plane  output  and  detection  of 
exploitable  cycle  frequency  features:  a)  .01  %  smoothing,  b)  5%  smoothing.  -lOdB 
SNR  input,  +2dB  in-band  SNR. 


Figure  5-7  Effect  ot  smoothing  on  bifrequency  plane  output  and  detection  of 
exploitable  cycle  frequency  features;  a)  .01  %  smoothing,  b)  5%  smoothing.  -15dB 
SNR  input,  -3dB  in-band  SNR. 
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Figure  5-8  Effect  of  smoothing  on  SPECCOA  output;  a)  0.005  % ,  b)  0.09% ,  c)  0.15%, 
d)  0.39%  with  resolution  products  of  1,  2,  4,  and  8,  respectively.  Input  SNR  was  -5dB 
with  a  chip  rate  of  0.0625  resulting  in  an  in-band  SNR  of  -I-  12dB. 
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Figure  5-9  Effect  of  decreasing  signal  to  noise  ratio  (SNR)  on  SPECCOA  output: 
a)  OdB,  b)  -SdB,  c)  -lOdB,  d)  -12dB.  Cycle  frequency,  0.0625,  carrier  frequency,  0.16, 
smoothing,  0.01%,  sample  lengths,  1920,  1920,  3840,  and  15,360,  respectively. 
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Figure  5-10  SPECCOA  output  results  with  various  interfering  signal  conditions,  a.,  = 
0.0625/T,. 
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VI.  CONCLUSIONS 


The  purpose  of  this  thesis  was  to  implement  the  SPECCOA  algorithm.  That 
objective  was  met.  Tests  conducted  showed  that  the  algorithm  was  successful  in 
producing  a  time-delay  output  in  simulated  environments  of  negative  signal  to  ratio 
(SNR)  and  interfering  signals.  These  tests  were  conducted  with  relatively  short  collect 
lengths  of  less  than  16,384  samples.  With  larger  collect  lengths  it  is  anticipated  that 
significantly  lower  SNR’s  can  be  achieved.  There  is  a  limit,  of  course,  as  to  how  long 
a  collect  length  can  be  gathered.  The  longer  the  collect  length,  the  longer  the  processing 
time.  For  systems  that  operate  in  near  real  time,  this  could  become  a  significant 
drawback.  Thus,  a  trade-off  must  be  made  between  the  time  necessary  to  process  the 
data  and  the  level  of  SNR  that  can  be  examined. 

The  most  significant  result  of  the  testing  was  that  a  high  degree  of  smoothing  was 
not  necessary  in  order  to  produce  a  time-delay  output.  In  fact,  a  high  degree  of 
smoothing  actually  degraded  the  performance  of  the  algorithm.  Very  small  values  of 
smoothing  were  required,  on  the  order  of  0.01  %  for  collect  lengths  of  less  than  10.(X)0. 
This  produced  a  resolution  product  of  near  unity.*  Significant  smoothing  was  still 
required,  however,  for  developing  plots  of  the  bifrequency  plane  in  order  to  see  what 
features,  if  any,  were  exploitable. 


'Recall  the  resolution  product.  AtAf.  equals  (%  smoothing/ 100)( sample  length)-r(  1-cn,). 
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Further  tests  are  needed  in  a  real-world  environment.  Data  exists  from  a  special 
experiment  that  was  conducted  in  a  real-world  environment  and  is  available  for 
processing  and  insertion  into  the  SPECCOA  algorithm. 

Finally,  other  methods  exist  to  calculate  the  spectral  correlation  densities  (SCDs) 
that  are  used  by  the  SPECCOA  algorithm.  Examples  are  the  FFT  Accumulation  Method 
(FAM)  and  the  Strip  Spectral  Correlation  Algorithm  (SSCA)  [Ref  11].  Testing  the 
SPECCOA  algorithm  using  other  methods  is  necessary  to  determine  if  any  modifications 
to  the  algorithm  are  required.  This  would  produce  a  general  purpose  SPECCOA 
algorithm  useful  with  a  variety  of  SCD  calculation  methods. 
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APPENDIX  A.  HYPERBOLA  DEVELOPMENT 


Following  [Refs.  2  and  9],  a  proof  is  given  that  the  geolocation  scenario  described 
in  Chapter  HI,  part  B  does  indeed  produce  a  hyperbola.  Figure  A-1  depicts  the  basic 
scenario.  Receiver  platforms  one  and  two  are  initially  located  along  the  x-axis  at 
coordinates  (-c,0)  and 
(c,0).  The  emitter,  located 
at  coordinates  (jr.y), 
transmits  a  spherical  wave 
that  travels  over  distances 
d,  and  dj  to  arrive  at  times 
t,  and  tj  at  receivers  one 
and  two,  respectively.  The 
time  difference  of  arrival, 
or  TDOA,  is  then 

TDOA  =  Ar  =  tj-r,  =  ^  -  (A.la) 

V  v 


Figure  A-1  Basic  TDOA  geometry. 


(A.lb) 
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where  v  is  the  propagation  speed  of  the  electromagnetic  wave.  Multiplying  both  sides 
by  V  to  get  an  equation  in  the  form  of  a  difference  of  distances  gives 

vAr  =  D  =  I  ^(jt  +c)^  I  .  (A. 2) 

Since  a  hyperbola  is  defined  as  the  set  of  all  points  whose  difference  in  distance  from  two 
fixed  points  is  constant,  manipulating  (A-2)  will  yield  the  standard  form  of  a  hyperbola. 
First,  moving  the  second  radical  to  the  left  side  of  the  equation,  squaring  both  sides,  and 
rearranging  gives 

+  (x-cf  +  *  2D\l{x-cf  +y’  =  (x*cf  +  y“, 

expanding  and  simplifying  leads  to 

4xc  -  =  lDsj(x-cfi  *  y^ .  (A.4) 

Again,  squaring  both  sides  and  rearranging  yields 

AD^c^  -  D*  =  16j:^c^  -  Ax^D^  -  AD^y" .  (A. 5) 

One  further  simplification  gives  the  final  result  of 

-4 - 

D^A  (4c ^  -  D^)/4 

which  is  the  standard  form  of  the  hyperbola: 
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'> 


(A.7) 


with  x-intercept  ±a  or  ±D/2.  and  asymptotic  to  the  lines  y  =  ± 


or 
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APPENDIX  B.  TEST  CASE  GRAPHS:  ao  EQUALS  CHIP  FEATURE 


TABLE  B1  PARAMETERS  FOR  GRAPHS  OF  HGURE  B-1:  VARYING  SNR 


input 

delay 

% 

smoothing 

resolution 

product' 

sample 

length 

chip 

feature" 

carrier 

frequency" 

SNR 

dB 

a) 

300 

0.01 

1 

1920 

0.0625 

0.16 

-0 

b) 

300 

0.01 

1 

1920 

0.0625 

0.16 

-5 

c) 

300 

0.01 

1 

3840 

0.0625 

0.16 

-10 

d) 

300 

0.001 

1 

15360 

0.0625 

0.16 

-12 

’AtAf  =  {%  smooth/ 100)(sample  length)/(l  -  aj  *  f,/f. 


TABLE  B2  PARAMETERS  FOR  GRAPHS  OF  HGURE  B-2:  VARYING  DELAY 


■ 

input 

delay 

% 

smoothing 

resolution 

product' 

sample 

length 

ofo,  chip 
feature" 

carrier 

frequency" 

SNR 

dB 

a) 

150 

0.01 

1 

1920 

0.0625 

0.16 

-5 

b) 

300 

0.01 

1 

1920 

0.0625 

0.16 

-5 

c) 

700 

0.01 

1 

1920 

0.0625 

0.16 

-5 

d) 

950 

O.OI 

1 

1920 

0.0625 

0.16 

-5 

'AtAf  =  (%  smooth/ 100)(sample  length)/(l  -  a„)  'f./f, 

TABLE  B3  PARAMETERS  FOR  GRAPHS  OF  HGURE  B-3: 
VARYING  SMOOTHING 


input 

delay 

% 

smoothing 

resolution 

product' 

sample 

length 

oTo,  chip 
feature" 

carrier 

frequency" 

SNR 

dB 

a) 

150 

0.005 

1 

1920 

0.0625 

0.16 

-5 

b) 

150 

0.09 

2 

1920 

0.0625 

0.16 

-5 

c) 

150 

0.15 

4 

1920 

0.0625 

0.16 

-5 

d) 

150 

0.39 

8 

1920 

0.0625 

0.16 

-5 

'AtAf  =  (%  smooth/ 100)(sample  length)/(l  -  a„)  '  f,/f. 
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TABLE  B4  PARAMETERS  FOR  GRAPHS  OF  HGURE  B-4:  VARYING  f 


input 

delay 

% 

smoothing 

resolution 

product' 

sample 

length 

do,  chip 
feature* 

carrier 

frequency* 

SNR 

dB 

a) 

150 

0.01 

1 

1920 

0.0625 

0.16 

-5 

b) 

150 

0.01 

1 

1945 

0.05 

0.16 

-5 

c) 

150 

0.01 

1 

1962 

0.04167 

0.16 

-5 

d) 

150 

0.01 

1 

1979 

0.03333 

0.16 

-5 

'AtAf  =  (%  smooth/ 100)(sample  length)/(l  -  aj  '  f,/f. 


TABLE  B5  PARAMETERS  FOR  GRAPHS  OF  HGURE  B-5:  VARING 


input 

delay 

% 

smoothing 

resolution 

product' 

sample 

length 

cto,  chip 
feature^ 

carrier 

frequency* 

SNR 

dB 

a) 

150 

0.01 

1 

1920 

0.0625 

0.10 

-5 

b) 

150 

0.01 

1 

1945 

0.0625 

0.16 

-5 

c) 

150 

0.01 

1 

1962 

0.0625 

0.30 

-5 

d) 

150 

0.01 

1 

1979 

0.0625 

0.42 

-5 

'AtAf  =  (%  smooth/ lOOX sample  length)/(l  -  aJ  -f,/f. 
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Figure  B-1  Effect  of  decreasing  signal  to  noise  ratio  (SNR)  on  SPECCOA  output: 
a)  OdB,  b)  -.5dB,  c)  -lOdB,  d)  -12dB.  Cycle  frequency,  .ObZS/T,,  carrier  frequency. 
.16/T„  smoothing,  0.01%. 
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Figure  B-2  Output  of  SPECCOA  algorithm  for  increasing  delays. 
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Figure  B-3  Effect  of  smoothing  on  SPECCOA  output:  a)  0.005%,  b)  0.09%.  c) 
0.15%,  d)  0.39%  with  resolution  products  of  1,  2,  4,  and  8,  respectively.  Carrier 
frequency  was  .16/T,. 
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Figure  B-4  Effect  on  SPECCOA  output  of  varying  the  chip  frequency,  f,.  Carrier 
frequency,  0. 16/T,. 
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APPENDIX  C.  TEST  CASE  GRAPHS:  ««  EQUALS  TWICE  CARRIER  FREQ. 


TABLE  Cl  PARAMETERS  FOR  GRAPHS  OF  HGURE  C-1:  VARYING  SNR 


■ 

input 

delay 

% 

smoothing 

resolution 

product’ 

sample 

length 

carrier 

frequency- 

SNR 

dB 

a) 

150 

0.20 

5 

1392 

0.32 

0.16 

0 

b) 

150 

0.20 

5 

1392 

0.32 

0.16 

-5 

c) 

150 

0.20 

5 

2392 

0.32 

0.16 

-10 

d) 

150 

0.0305 

5 

11,141 

0.32 

0.16 

-15 

'AtAf  =  (%  smooth/ 100)( sample  length)/(l  -  a„)  '  f,/f„  =  0.0625/Ts. 


TABLE  C2  PARAMETERS  FOR  GRAPHS  OF  HGURE  C-2:  VARYING  DELAY 


■ 

input 

delay 

% 

smoothing 

resolution 

product' 

sample 

length 

carrier 

frequency* 

SNR 

dB 

a) 

150 

0.20 

5 

1392 

0.32 

0.16 

-5 

b) 

300 

0.20 

5 

2785 

0.32 

0.16 

-5 

c) 

700 

0.05 

5 

5570 

0.32 

0.16 

-5 

d) 

900 

0.05 

5 

5570 

0.32 

0.16 

-5 

’AtAf  =  (%  smooth/ 100)(sample  Iength)/(1  -  a„)  *  f,/f„  =  O.O625/T5. 
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TABLE  C3  PARAMETERS  FOR  GRAPHS  OF  HGURE  C-3: 
VARYING  SMOOTHING 


input 

% 

resolution 

sample 

carrier 

SNR 

delay 

smoothing 

product' 

length 

frequency- 

dB 

a) 

150 

0.01 

1 

1392 

0.32 

0.16 

-5 

b) 

150 

0.195 

4 

1392 

0.32 

0.16 

-5 

c) 

150 

0.35 

8 

1392 

0.32 

0.16 

-5 

d) 

150 

0.50 

11 

1392 

0.32 

0.16 

-5 

*AtAf  =  (%  smooth/ 100)(sample  length)/(l  -  «„)  '  f,/f„  =  O.O625/T5. 


TABLE  C4  PARAMETERS  FOR  GRAPHS  OF  HGURE  C-4:  VARYING 


■ 

input 

delay 

% 

smoothing 

resolution 

product' 

sample 

length 

carrier 

frequency- 

SNR 

dB 

a) 

150 

0.20 

5 

1638 

0.20 

0.10 

-5 

b) 

150 

0.20 

5 

1515 

0.21 

0.13 

-5 

c) 

150 

0.20 

5 

1392 

0.32 

0.16 

-5 

d) 

150 

0.20 

5 

1269 

0.38 

0.19 

-5 

'AtAf  =  (%  smooth/ I00)(sample  length)/(I  -  a„)  '  f,/f,,  =  0,0625/T,. 
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Figure  C-1  Effect  of  decreasing  signal  to  noise  ratio  (SNR)  on  SPECCOA  output: 
a)  OdB,  b)  -5dB,  c)  -lOdB,  d)  -I5dB.  Cycle  frequency,  a  =  carrier 
frequency,  .16/T„  AtAf  =  5. 
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APPENDIX  D.  TEST  CASE  GRAPHS:  EQUALS  TWICE  CARRIER  FREQ. 

PLUS  CHIP  FEATLTIE 


TABLE  D1  PARAMETERS  FOR  GRAPHS  OF  HGURE  D-1:  VARYING  SNR 


■ 

input 

delay 

% 

smoothing 

resolution 

product* 

sample 

length 

carrier 

frequency* 

SNR 

dB 

a) 

150 

0.02 

1 

1264 

0.3825 

0.16 

0 

b) 

150 

0.02 

1 

1264 

0.3825 

0.16 

-5 

c) 

150 

0.02 

1 

2529 

0.3825 

0.16 

-10 

d) 

150 

0.005 

1 

10,117 

0.3825 

0.16 

-15 

’iAtAf  =  {%  smooth/ 100)(sample  length)/(l  -  o(„)  ’  f,/T„  =  0.0625/T,. 


TABLE  D2  PARAMETERS  FOR  GRAPHS  OF  HGURE  D-2:  VARYING  DELAY 


■ 

input 

delay 

% 

smoothing 

resolution 

product* 

sample 

length 

carrier 

frequency* 

SNR 

dB 

a) 

150 

0.02 

1 

2529 

0.3825 

0.16 

-5 

b) 

300 

0.02 

1 

2529 

0.3825 

0.16 

-5 

c) 

700 

0.02 

1 

2529 

0.3825 

0.16 

-5 

d) 

900 

0.02 

1 

2529 

0.3825 

0.16 

-5 

'AtAf  =  (%  smooth/ 100)(sample  length)/(l  -  a„)  ’  f,/T,.  =  0.0625/T, 
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TABLE  D3  PARAMETERS  FOR  GRAPHS  OF  HGURE  D-3; 
VARYING  SMOOTHING 


■ 

input 

delay 

% 

smoothing 

resolution 

product' 

sample 

length 

carrier 

frequency* 

SNR 

dB 

a) 

150 

0.05 

2 

1264 

0.3825 

0.16 

-5 

b) 

150 

0.10 

3 

1264 

0.3825 

0.16 

-5 

c) 

150 

0.30 

7 

1264 

0.3825 

0.16 

-5 

d) 

150 

0.50 

11 

1264 

0.3825 

0.16 

-5 

'AtAf  =  (%  smooth/ 100)(sample  length)/(l  -  aj  '  f,/T,.  =  0.0625/T.. 


TABLE  D4  PARAMETERS  FOR  GRAPHS  OF  HGURE  D-4:  VARYING  f. 


■ 

input 

delay 

% 

smoothing 

resolution 

product* 

sample 

length 

“0^ 

carrier 

frequency* 

SNR 

dB 

a) 

150 

0.02 

1 

1510 

0.2625 

0.10 

-5 

b) 

150 

0.02 

1 

1387 

0.3225 

0.13 

-5 

c) 

150 

0.02 

1 

1264 

0.3825 

0.16 

-5 

d) 

150 

0.02 

1 

1141 

0.4425 

0.19 

-5 

'AtAf  =  (%  smooth/ lOOKsample  Iength)/(1  -  a„)  -  f,/T„  f,  =  O.O625/T5. 
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magnitude  maqnili 


-500  —4.00  —300  —200  -  i  00  0 

tou  /  Ts 

ae<ay...l50  sannoiea  output  paok . i  50  so'^pies 


(a)  1264  samples 


Input  dalay,..l50  aampiaa  output  peak . 150  sompies 

(b)  1264  samples 


input  delay.  ..150  samples  output  oeok . i  50  sompies 


(c)  2529  samples 


input  delay...  150  samples  output  peak . 150  samples 

d)  10,1 17  samples 

Figure  D-1  Effect  of  decreasing  SNR  on  SPECCOA  output:  a)  OdB,  b)  -5dB.  c)  -lOdB. 
d)  -15dB.  Carrier  frequency,  0. 16/T„  f^,  0.0625/Ts,  =  2  f„  +  f„. 
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moqnitude  magnitude  magnitude  magnitude 


C.2h 


-1000  -800  -600  -4-00  -200 

tou  /  ■'« 

•  '^out  cl*iay...i50  samDi*s  cutout  oaok . 'SO  sc^o'** 


input  C«loy...300  acrrioiaa  cutout  oacK . 300  acnnpiaa 


-1000 


-800 


inout  daiay...1’00  acmolas 


-600  -AOO  -200 

tou/Ts 

cutout  Deck . ~C0  san-iQias 


Figure  D-2  Output  of  SPECCOA  algorithm  for  increasing  delays. 
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S' 

£ 


OJ 

■o 

3 


50 


-500 


input  deioy...l50  samples 


output  peck . -^7  sompies 


Figure  D-3  Effect  of  smoothing  on  SPECCOA  output:  a)  0.05%,  b)  0.10%.  c)  0.30%. 
d)  0.50%.  with  resolution  products  of  2,  3,  7,  11,  respectively.  Carrier  frequency, 
0.1 6/T,. 
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-500  -4.00  -300  -200  -^00 

tau/ Ts 


'  no  uX.  a«iay...‘^50  somoiea  outDut  D«ok . tSO  samoias 

(a)f„,  0.1/T,  0.2625/T,. 


input  aeiay...150  sompias  output  ceok . 150  sompies 


(b)f„0.13/T,.  a„  0.3225/T, 


(c)  fo,  0.I6/T,.  a„,  0.3825/T,. 


(d)f„,0.1/T,.  a„,  0.4425/T,. 

Figure  D-4  Effect  of  varying  the  carrier  frequency,  f^,  on  SPECCOA  output. 
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APPENDIX  E.  COMPUTER  PROGRAMS 


This  appendix  contains  four  computer  programs  and  routines  that  were  written 
during  the  course  of  working  on  this  paper:  tdoasmth.bat,  speca.c,  fcnvn.c,  and 
plotdata.m.  The  main  algorithm  written  is  called  speca  and  has  a  companion  program, 
fcnvrt.  that  is  used  to  adjust  the  data  file  headers  obtained  from  the  SSPI  software 
package  so  that  they  can  be  used  by  speca.  If  another  method  is  used  to  calculate  the 
spectral  correlation  density  functions  needed  in  the  SPECCOA  algorithm,  then  their 
output  data  files  need  to  be  formatted  properly  as  noted  in  the  speca  program,  tdoasmth 
is  a  controlling  script  file  that  calls  the  signal  generating  program,  spectral  correlation 
density  calculating  program,  and  the  TDOA  calculating  algorithm,  plotdata  takes  the 
output  data  produced  by  the  speca  program  and  plots  it.  Each  program  is  thought  to  be 
well  commented  so  that  an  additional  in  depth  discussion  will  not  be  given  here.  Two 
areas  requiring  some  comment,  however,  are  generating  the  test  signals  and  the  spectral 
correlation  densities,  both  of  which  use  the  SSPI  software  package. 

A.  SIGNAL  GENERATION  SOFTWARE 

This  section  describes  the  software  created  by  Statistical  Signal  Processing,  Inc. 
that  was  used  to  simulate  signals  to  test  the  SPECCOA  algorithm.  Table  E-1  lists  the 
parameters  that  were  varied  to  obtain  the  desired  characteristics  of  the  test  signal.  One 
or  more  of  the  parameters  shown  in  bold  was  varied  for  each  test  run.  The  other 
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parameters  could  also  be  varied  but  were  essentially  left  in  their  default  setting  when  the 
input  file  was  first  generated.  More  information  can  be  obtained  from  [Ref.  10]  which 
should  be  consulted  frequently  when  using  the  SSPI  software. 


B.  GENERATING  SPECTRAL  CORRELATION  DENSITIES 

The  spectral  correlation  and  cross  spectral  correlation  densities  were  generated 
using  the  SSPI  software.  Script  file  tdoasmth.bat  was  used  to  call  the  SSPI  software 
package  to  perform  the  calculations.  The  steps  involved  for  a  particular  test  run  were 
as  follows: 

•  adjust  parameters  in  the  file  pam Jilnam.inp  to  create  the  desired  test  signal. 

•  set  Qfo  in  the  alphas  file  to  the  value  of  cycle  frequency  being  exploited. 

•  adjust  smoothing  and  cycle  frequency  as  necessary  in  the  tdoasmth.bat  file; 
smoothing  is  command  input  for  the  sxaf  command  (e.g.,  -w  0.01),  and  the  alpha 
value  is  a  command  line  input  for  the  speca  routine. 

•  call  tdoasmth.bat  to  calculate  the  spectral  and  cross-spectral  densities;  output  files 
are  surjyxj)ut  and  surfx  out. 

•  generated  data  files  can  now  be  processed  by  the  SPECCOA  algorithm  to 
determine  the  time  delay;  algorithm  output  is  stored  in  Syxx.dat  file  which  was 
read  directly  into  the  MATLAB  routine plotdata.m  to  plot  the  SPECCOA  algorithm 
output. 
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TABLE  E-1  INPUT  PARAMETERS  IN 
SIGNAL  GENERATION  PROGRAM 


nsamples 

2048 

samples  per  symbol 

16 

bits  per  symbol 

1 

constellation 

psk 

pulse  type 

nyq 

pulse  BW 

l.OOOOOOe+00 

pulse  tail  exp 

-2 

delay  samples 

300 

real  freq  shift 

1 

carrier  freq 

0.16e+00 

carrier  phase  deg 

O.OOOOOOe+00 

over  sample 

1 

under  sample 

1 

signal  power  dB 

0.5000e+01 

noise  power  dB 

1.5000e+01 

ASCn  or  binary 

ASCn 

real  output 

0 

time  output 

1 

freq  output 

0 
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TDOASMTH.BAT 


#This  file  IS  a  batch/script  file  to  run  the  separate  programs  needed  to 
#generate  a  test  signal,  calculate  the  spectral  and  cross-spectral  correlation 
#densities,  and  determine  the  time  difference  of  arrival  using  the  SPECCOA 
#algorithm.  Refer  to  the  SSPI  software  manual  for  details  of  'pam',  'combine', 
#and  ' sxaf '  programs. 

#create  the  test  signals  x  and  y.  y  =  x(t-delay) 

#the  output  files  generated  are  pam_xsig.tim,  pam_ysig.tim 
"^modifying  the  pam_xsig.inp  or  pam_ysig.inp  file  for  binary  output 
#vice  ASCII  output  will  make  the  program  run  faster 

pam  -bits  -1534  xsigo  #original  signal  before  noise  added 

pam  -bits  -1534  ysigo  #original  signal  before  noise  added 

pam  -noise  -2066  nnoise  #generate  noise  with  independent  seed 

pam  -noise  -4065  mnoise  #generate  m-noise  with  different  ind  seed 

combine  pam_xsigo . tim  pam_nnoise .tim  pam_xsig.tim  #add  nnoise  to  xsig 
combine  pam_ysigo ,tim  pam_mnoise .tim  pam_ysig.tim  #add  mnoise  to  ysig 

tcross  correlate  xsig,  ysig  to  generate  the  cross  spectral  correlation,  Syx. 
fusing  the  (-f  a)  option  makes  the  output  file  ASCII,  (-f  b)  for  binary, 
fuse  ' fcnvrt'  program  to  make  the  output  file  header  compatible  with 
fthe  ' speca'  program.  NOTE:  cycle  frequency,  alpha,  is  read  from  the  default 
ffile  named  'alphas'. 

sxaf  -f  a  -w  0.01  -o  surfyx_out  pam_ysig.tim  pam_xsig.tim 
fcnvrt  surfyx_out  syx_out  f  fcnvrt  infile  outfile 

fautocorrelate  xsig  to  generate  Sx.  'fcnvrt'  makes  output  file  header 
f compatible  with  the  'speca'  program. 

sxaf  -f  a  -w  0.01  -o  surfx_out  pam_xsig.tira 
fcnvrt  surfx_out  sx_out 

fmultiply  (Syx)  &  con jugate (Sx) .  Inverse  Fourier  transform  the  product. 
fThe  specific  cycle  frequency,  alpha,  is  an  input  parameter;  0.0625  is  shown. 
fThe  result  is  stored  in  the  data  file  'Syxx.dat',  which  can  be  loaded  directly 
finto  MATLAB  for  plotting.  Format  is  a  col  of  real  nos.  and  a  col  of  imag.  nos. 

speca  0.0625  syx_out  sx_out  Syxx.dat 
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SPECA.C 


#include  <stdlib.h> 

#include  <stdio.h> 

#include  <math.h> 

#include  "/tools2/SSPI/pam/f ft . c" 

#include  " /tools2/SSPI/pam/radix . c" 

main(argc,  argv) 

int  argc; 
char  *argv [ ] ; 

{ 

int  1,  r_]_flagl,  num_data_ptsl,  n,  n2,  N,  pwr2; 

int  r_]_flag2,  r_]_f lag_out,  num_data_pts2,  nuin_of_loops,  zeroref,  tau; 

char  *file2,  *wfile; 

float  alpha,  rtempl,  rtemp2,  itempl,  itemp2; 

COMPLEX  *Syxx; 
static  COMPLEX  sigtemp; 

FILE  *fprl,*fpr2,  *fpw; 

1=0; 

alpha  =  atoi (argv [ 1 ] ) ;  /*  input  cycle  frequency  parameter  to  used  in  SPECCOA  algorithm 

filel  =  argv [2];  /*  Syx  data  file 

file2  =  argv [3];  /*  Sx  data  file 

wfile  =  argv[41;  I*  Syxx  output  data  file 


/*  »/ 

/*  This  program  computes  the  TIME  DIFFERENCE  OF  ARRIVAL  (TDOA)  from  */ 

/*  two  input  files,  one  a  delayed  version  of  the  other  */ 

/*  */ 

/*  written  by  LT  Timothy  A.  Benson  ♦/ 

/*  December  1992  */ 

/*  ♦/ 

/* - */ 

/»  files  have  the  following  format 


input  files  - 

real  or  complex  flag 
num  data  points 
data : 

real_part  imag_part 

output  file  - 

real  or  complex  flag 
num  data  points 
data : 

real_part  imag_part 

-  */ 

/*  */ 
/♦*♦**■••♦♦**♦***♦*♦♦♦»*♦♦♦♦****»♦*♦*♦******♦,******♦♦«*******♦♦**♦♦*, ♦*«****^ 
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VARIABLE  DEFINITIONS 


*  '* 


*  / 
♦/ 
*  / 


/  * 
/* 


*/ 


Syxx 

1 

r_3_f lag 
num_data_pts 
n 
N 

pwr2 

n2 

num_lQops 

zeroref 
tau 
alpha 
rtempl, 2 
itempl, 2 
sigtemp 


-  contains  the  final  data  of  this  program 

-  a  counter 

-  real-imaginary  flag  for  files  1,  2,  and  output  file 

-  number  of  data  pts  for  file  1,  2 

-  equals  number  of  data  pts  in  smaller  file 

-  power  of  two  next  avbove  n 

-  a  number  such  that  10^pwr2  is  greater  than  or  equal  to  n 

-  n2  =  n/2 

-  loop  through  this  many  times  to  read  amount  of  data  from  both 
equal  to  the  amount  of  data  in  the  smaller  file 

-  finds  zero  index  for  data  length  -N/2  to  N/2 

-  any  value  between  -N/2  to  N/2.  Corresponds  to  the  lag  value 

-  cycle  frequency  of  spectral  correlation  being  processed 

-  real  varialbe  into  which  data  is  read  from  in-file 

-  imaginary  varialbe  into  which  data  is  read  from  in-file 

-  a  temporary  variable  to  hold  values  of  Syxx  during  manipulati 


/*  check  for  the  correct  no.  of  command  line  arguments  */ 
if  (argc  ==  1)  { 

printf ("speca  alpha  filel  file2  outfile") 
printfC  -  takes  spectral  data  from  filel  and  file2\n") 
printf ("and  computes  TDOA  using  SPECCOA  algorithm: \n") 
printf ("delay  =  real [ifft  (SyxSx*)exp( j*pi*alpha*tau) \n") ; 
exit  ( )  ; 

else 

if  (argc  !=  5)  { 

printf ( "speca  ...  improper  no.  of  arguments\n" ) ; 

printf ("proper  format  is:  speca  alpha  Syxfil  Sxfil  outfileXn"); 

printf  ( " \n" ) ; 

exit  ()  ; 


/*  open  files  and  assign  pointers  to  them  */ 
fprl  =  fopen( filel,  "r"); 
fpr2  =  fopen(file2,  "r"); 
fpw  =  f open (wf lie,  "w"); 


/*  determine  whether  data  file  has  complex  numbers  and  get  number  of  data  points  */ 
/*  first  file  */ 

fscanf (fprl, "%d  %d",  &r_3_flagl,  &num_data_ptsl) ; 

/*  second  file  */ 

fscanf (fpr2, "%d  %d",  &r_]_flag2,  &num_data_pts2) ; 
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*  loop  through  num  of  data  points  for  the  file  with  the  lesser  amount  */ 

if  (num_data_ptsl  <  num_data_pts2) 
n  =  num_data_ptsl; 

else 

n  =  num_data_pts2  ; 
nura_of_loops  =  n; 


/*  determine  if  n  is 
pwr2=radix (n) ;  /* 

N=l; 

for  (i=0;  i<  pwr2; 
N=2*N; 


a  power  of  two;  if  not  either  truncate  or  zero-pad  to  a  power 
set  equal  to  smallest  power  of  two  >=  n  */ 

++i) 


/*  determine  whether  to  zero-pad  or  truncate,  if  n  is  greater  than  35%  of  the  way 
to  the  next  power  of  2,  then  zero-pad,  otherwise  truncate  */ 

if(n  >  1.35*N/2)  { 

printf ("n=%d,  N=%d,  n  will  be  zero-padded  to  %d\n",  n,  N,  N)  ; 


if(n  <=  1.35*N/2)  { 

N=N/2; 

printf ("n=%d,  N=%d,  truncating  n  to  %d\n",  n,  N,  N) ; 


/*  test  line  */ 

printf ("1;  n=%d,  N=%d\n",  n,  N) ; 

/*  allocate  space  for  Syxx  vector  */ 

Syxx= (COMPLEX*) calloc (N,  sizeof (COMPLEX) )  ; 
if  (Syxx  ==  NULL)  { 

printf ("speca ...  insufficient  space  to  allocate  Syxx\n"); 
exit  ()  ; 

) 


/*  if  n  IS  not  a  power  of  2,  then  zero  pad  symetrically  up  to  a  power  of  2,  N  */ 
if  (N  >  n)  { 

/*  determine  the  number  of  zeros  that  have  to  be  padded  */ 
numzer  =  N  -  n; 

) 

/*  test  line  */ 

printf ("2:  n=%d,  N=%d\n",  n,  N) ; 
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/ 


OOP  to  read  in  all  data  and  conjugate  and  multiply  as  necessary  to 
orm  product  of  (Syx)  *  CQn]ugate (Sx) 

for (1=0;  i<N;  ++i)  { 

/*  if  zero-padding  needed  (le  N  >  n)  */ 
if(N  >  n)  i 

/*  zero-pad  the  first  (num.  of  zeros) /2  elements  of  Syxx  */ 
if  (i  <  (numzer+1 ) /2 )  { 

Syxx [i] .r  =  0; 

Syxx [i] . 1  =  0; 


/*  zero-pad  the  last  (num.  of  zeros) /2  elements  of  Syxx  */ 
if  (i  >=  (N  -  numzer/2))  { 

Syxx [ 1 ] . r  =  0; 

Syxx [i] . 1  =  0; 


}  /*  end  loop  for  zero-padding  '/ 


if  (r_]_flagl  ==  1  &&  r_j_flag2  ==1)  {  /*  both  files  are  REAL  only 

r_]_flag_out  =  1; 
fscanf (fprl, "%e" ,  irtempl) ; 
fscanf (fpr2, "%e",  &rtemp2); 

Syxx[i].r  =  (rtempl)  *  (rtemp2); 

Syxx [i] . 1  =  0; 

} 

if  (r_j_flagl  ==  1  &&  r_]_flag2  ==2)  {  /*  filel  real,  file2  complx  * 

r_3_flag_out  =  2; 
fscanf (fprl, "%e",  &rtempl); 
fscanf ( fpr2, " %e%e" ,  &rtemp2,  &itemp2); 

Syxx[i].r  =  rtempl  *  rtemp2; 

Syxx[i].i  =  rtempl  *  (-1.0)  *  itemp2; 

} 


if  (r_]_flagl  ==  2  &&  r_j_flag2  ==1)  {  /*  filel  complx,  file2  real  */ 

r_j_flag_out  =  2; 

fscanf (fprl, "%e%e",  &rtempl,  iitempl); 
fscanf (fpr2,  "%e",  &rtemp2); 

Syxx[i].r  =  rtempl  *  rtemp2; 

Syxx[i].i  =  rtemp2  *  itempl; 
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/*  add  data  to  Syxx  if  i  is  between  zero-padded  ends 

if  (i  >=  (numzer  +l)/2  &&  i  <  (N  -  numzer/2))  { 

if  (r_]_flagl  ==  2  &&  r_]_flag2  ==2)  {  /*  both  files  are  complx  */ 

r_3_flag_out  =  2; 

f scanf ( fprl, " %e%e" ,  srtempl,  iitempl); 
fscanf (fpr2, "%e%e",  &rtemp2,  &itemp2); 

Syxx[i].r  =  (rtempl  *  rtemp2  -  itempl  ♦  (-1.0)  ’  itemp2)/n; 
Syxx[i].i  =  (rtempl  *  (-1.0)  *  itemp2  +  rtemp2  *  itempD/n; 

} 

} 

)  /*  end  of  for  loop.  all  data  is  multiplied  and  stored  in  Syxx  */ 


/*•  test  line  */ 

printfC'all  data  has  been  read  and  multiplied\n" ) ; 


/*  test  line  */ 

printf ("beginning  fft  operationXn" ) ; 

/*  inverse  Fourier  transform  Syxx  */ 

pwr2=radix (n) ;  /*  set  equal  to  smallest  power  of  two  >=  n  */ 

N=l; 

for  (i=0;  i<  pwr2;  ++i) 

N=2*N; 


/*  test  line  */ 
printf ("n=%d,  N=%d\n",  n,  N) ; 

fft (Syxx,  n,  -1,  2); 

/*  Syxx  -  COMPLEX  array  to  be  transformed 

n  -  number  of  points  in  signal  -  power  of  2  >=  n 

-1  -  inverse  FFT 

2  -  normalized  FFT 


/*  shift  result  from  (0  to  N-1)  to  (-N/2  to  N/2)  */ 
for  (i=0;  i  <  n2;  ++i)  { 

sigtemp  =  *(Syxx  +  i); 

*(Syxx  +  i)  =  *(Syxx  +  i  +  n2)  ; 

*(Syxx  +  1  +  n2)  =  sigtemp; 
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/*  test  line  */ 

printfC'ifft  routine  finished\n") 
printf  ("multiplying  by  final  factorin’’) 

/*  multiply  by  exp ( ] *pi*alpha*tau)  where  zero  ref  is  n/2  +  1  */ 

/*  this  expression  is  equivalent  to  cos (pi*alpha’tau)  +  ] *sin  (pi’alpha'tau)  */ 
zeroref  =  n/2  +  1; 
for  (i=0;  i<n;  ++i)  { 

tau  =  -(zeroref  -  i) ; 

Syxx[i].r  =  Syxx[i].r  *  cos  (PI*alpha*tau)  -  Syxx[i].i  *  sin  (PI *alpha*tau) ; 
Syxx[i].i  =  Syxx[i].r  *  sin  (PI*alpha*tau)  +  Syxx[i].i  ♦  cos  (PI*alpha’'tau)  ; 

} 


/*  print  data  to  file  */ 

printf  ("writing  data  to  filein’’); 
for  (i=0;  i<n;  ++i) 

fprintf  ( fpw,  "  %e  %e\n’’,  Syxx[i].r,  Syxx[i].i); 
printf  ( "data  is  written  to  file\n’’); 


/♦  close  files  ♦/ 
fclose (fprl) ; 
f close (fpr2) ; 
fclose ( fpw)  ; 


}  /*  end  of  program  */ 
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FCNVRT . C 


>finclude  <stcilib.h> 

=Jinclude  <stdio.h> 

main (argc,  argv) 

int  argc; 
char  *argv[]; 

{ 

int  1,  r_3_flag,  num_data_pts,  tempint; 

char  *rfile,  *wfile; 

float  real,  imag,  tempfloat; 

FILE  *fpr,  ’fpw; 
rf lie  =  argv [ 1 ] ; 
wfile  =  argv [2]; 


/*  */ 

/♦  This  program  converts  a  data  file's  header  to  run  with  ' speca'  ♦/ 

/*  */ 

/*  written  by  Timothy  A.  Benson  */ 

/*  December  1992  */ 

/*  ’/ 

/***♦♦**»*•***»*♦****♦»***♦**♦******»*******»***♦*»*»*****♦»*♦*<♦*»♦♦*♦*♦/ 

/* - ./ 

/*  ’/ 

/*  file  to  be  read  has  the  following  format  -  note  line  spacing. 

(see  page  sxaf-8  in  SSPI  manual) 


<...>  -  indicates  nos.  read  by  program  to  control  program  flow, 
i.e.,  other  nos.  are  not  needed  but  correct  linespacing  is  required. 

<real  or  complex  flag> 

<num  cuts>  alpha  min  alpha  max 

num  data  points  max  fmin  max  fmax  max 

<alpha>  <numf>  min  fmax 

data 

-  ♦/ 
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/♦*»»»******»*»********»*********************»»»»**»*****»»**»»*»»-»*»***y' 
/*  */ 
/♦  VARIABLE  DEFINITIONS  ♦/ 

1  -  a  counter 

r_]_flag  -  real-  imaginary  flag 

num_dat_pts  -  number  of  data  points  in  in-file 

tempint  -  a  temporary  variable 

tempfloat  -  a  temporary  variable 

real  -  data  variables 

imag  -  data  variables 


/*  check  for  the  correct  no.  of  command  line  arguments  */ 
if  (argc  ==  1)  i 

printf ( " f cnvrt  infile  outfile  -  converts  header  format  for  use  by  fmultXn"); 
exit  0  ; 

} 

else 

if  (argc  '.  =  3)  { 

printf ("fcnvrt  ...  improper  no.  of  arguments \n" ) ; 
printf ("proper  format  is:  fcnvrt  infile  outfileVn"); 
printf ("\n") ; 
exit  0  ; 


/*  open  files  and  assign  pointers  to  them  */ 
fpr  =  fopen(rfile,  "r"); 
fpw  =  fopen(wfile,  "w"); 


/*  determine  whether  data  file  has  complex  numbers  and  the 
total  number  of  data  points.  */ 

fscanf (fpr,  "%d%d%e%e",  &r_]_flag,  itempint,  &tempfloat,  &tempfloat); 
fscanf (fpr, "%d%e%e",  &tempint,  stempfloat,  &tempf loat) ; 
fscanf (fpr, "%e%d",  &tempfloat,  &num_datajpts) ; 


/*  print  header  in  output  file  for  data  information  */ 
fprintf (fpw, "%d\n%d\n",  r_]_flag,  num_data_pts) ; 


/*  adjust  file  pointer  to  first  line  of  data  */ 
fscanf  (fpr, "%e%e",  &tempfloat,  stempfloat); 
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o 


'  ■*  loop  oo  copy  data  into  outfile  *  ■' 


1  =  0; 

while  (1  <  num_data_pts )  { 

+  +  1  ; 

if  (r_]_flag  ==  1)  | 

f scant ( fpr, " %e" ,  &real); 
fprintf (fpw, "%e\n",  real); 


else  if  (r_]_flag  ==  2)  { 

fscanf (fpr, "%e%e",  &real,  &imag) ; 
fprintf (fpw, "%e  le\n",  real. 


}  /*  end  while  loop  to  copy  all  data  */ 


lose  files  */ 
close ( fpr) ; 
f close ( fpw) ; 


imag)  ; 
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PLOTDATA.M 


%This  program  plots  data  from  a  file  named  Syxx.dat  created  by  tdoa.bat 


clear 

clg 

delay=101;  %value  of  delay  that  was  inserted  into  test  signal 

load  Syxx.dat 
len=length (Syxx) ; 

Syxxreal=  (Syxx  ( : , 1) ' ) ; 

Idetermine  lag  value  of  peak 
[maxval  indx] =max (Syxxreal)  ; 
lag=indx  -  len/2  -  1; 

%plot  data 

xaxis=-len/2 : len/2-1; 

V='[-len/2  0  -1 .2 *max  (-Syxxreal)  1 .2*max  (Syxxreal)  ]  ; 

V=[-512  0  -1 .2 *max (-Syxxreal)  1 .2*max (Syxxreal)]; 


axis (V) ; 
subplot  (211) 

plot (xaxis, Syxxreal) , grid 
title (' SPECCOA  OUTPUT') 
xlabel (' tau/Ts' ) 
ylabel ( '  magnitude' ) 

axis; 
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